// convex_hull.cpp
// 04/16/2011: created

#include <fstream>
#include <cassert>
#include <limits>
#include <algorithm>
#include <QtGui>
#include <QtCore/QtAlgorithms>
#include "geometric_problem.h"
#include "Log.h"

void ConvexHull::loadData(const std::string& fName)
{
  m_point.clear();
  m_vertexIdx.clear();

  int n, d;
  std::ifstream is(fName.c_str());
  is >> n >> d;

  Log::m_log << "n=" << n << " d=" << d << endl;
  for (int i = 0; i < n && !is.eof(); i++) {
    m_point.push_back(vec(d));
    vec& v = m_point.back();
    v.load(is, arma::arma_ascii);
    assert(v.n_elem == d);
    //Log::m_log << "v[" << i << "]" << endl;
  }
  is.close();
}

struct Compare {
  const std::vector<double> &m_cos_theta, m_dist;
  Compare(const std::vector<double>& cos_theta, const std::vector<double>& dist) : m_cos_theta(cos_theta), m_dist(dist) {}
  bool operator()(int i, int j) const {
    return m_cos_theta[i] < m_cos_theta[j] || ( fabs(m_cos_theta[i] - m_cos_theta[j]) < 1e-16 && m_dist[i] < m_dist[j] );
  }
};

// O(n log n) algorithm
void ConvexHull::solve()
{
  int n = m_point.size();

  // find left most point (if tight, use upper most)
  int pIdx = 0;
  for (int i = 1; i < n; i++) {
    vec& p = m_point[i];
    if (p[0] < m_point[pIdx][0] || (p[0] == m_point[pIdx][0] && p[1] < m_point[pIdx][1])) 
      pIdx = i;
  }

  // compute cos of angles between Oy and (X[i]-X[pIdx]) and then sort them
  std::vector<double> cos_theta(n);
  std::vector<double> dist(n);
  QVector<int> idx(n);
  for (int i = 0; i < n; i++) {
    vec v = m_point[i] - m_point[pIdx];
    double d = arma::norm(v, 2);
    if (d < 1e-16)
      cos_theta[i] = std::numeric_limits<double>::infinity();
    else
      cos_theta[i] = v[1] / d;
    dist[i] = d;
    idx[i] = i;
  }
  Log::m_log << "sorting ... "; Log::m_log.flush();
  //qSort(idx.begin(), idx.end(), Compare(cos_theta, dist));
  ::qSort(idx.begin(), idx.end(), Compare(cos_theta, dist));
  Log::m_log << "done" << endl;

  // main loop
  m_vertexIdx.clear();
  m_vertexIdx.push_back(pIdx);
  m_vertexIdx.push_back(idx[0]);
  int next = 1, current = 1;
  while (next < n) {
    if (cos_theta[idx[next]] == std::numeric_limits<double>::infinity()) break;
    int prev = current - 1;
    const vec& prevPoint = m_point[m_vertexIdx[prev]];
    const vec& currentPoint = m_point[m_vertexIdx[current]];
    const vec& nextPoint = m_point[idx[next]];
    vec v1 = currentPoint - prevPoint, v2 = nextPoint-currentPoint;
    if (v1.cross(v2) > 0) {
      m_vertexIdx.push_back(idx[next]);
      current++;
      next++;
      //Log::m_log << "add vertex: " << idx[next] << " (" << m_point[idx[next]][0]  << ", " << m_point[idx[next]][1] << ")" << endl;
    }
    else {
      //Log::m_log << "remove vertex: " << m_vertexIdx[current] << " (" << m_point[m_vertexIdx[current]][0]  << ", " << m_point[m_vertexIdx[current]][1] << ")" << endl;
      m_vertexIdx.pop_back();
      current--;
    }
  }
}

// O(n^2) algorithm
//void ConvexHull::solve()
//{
//  int n = m_point.size();
//  // find left most point (if tight, use upper most)
//  int pIdx = 0;
//  for (int i = 1; i < n; i++) {
//    vec& p = m_point[i];
//    if (p[0] < m_point[pIdx][0] || (p[0] == m_point[pIdx][0] && p[1] < m_point[pIdx][1])) 
//      pIdx = i;
//  }
//
//  std::vector<bool> used(n, false);
//  int start = pIdx;
//  do {
//    Log::m_log << "current vertex: " << pIdx << " (" << m_point[pIdx][0]  << ", " << m_point[pIdx][1] << ")" << endl;
//    m_vertexIdx.push_back(pIdx);
//    int nextIdx = -1;
//    double dist = std::numeric_limits<double>::infinity();
//    for (int i = 0; i < n; i++) {
//      if (i == pIdx || used[i]) continue;  // skip current point and used points
//      if (nextIdx == -1) nextIdx = i;
//      vec v1 = m_point[i] - m_point[pIdx], v2 = m_point[nextIdx]-m_point[pIdx];
//      double d = arma::norm(v1, 2);
//      double cross = v1.cross(v2);
//      if (cross < 0 || (cross < 1e-16 && dist > d)) {
//        nextIdx = i;
//        dist = d;
//      }
//    }
//    pIdx = nextIdx;
//    used[pIdx] = true;
//  } while (pIdx != -1 && pIdx != start);
//}

void ConvexHull::draw(Plotter& pl)
{
  QPen pen(Qt::SolidLine);
  pen.setWidth(2);

  pen.setColor(Qt::black);
  pl.setPen(pen);
  int n = m_point.size();
  for (int i = 0; i < n; i++) {
    vec& p = m_point[i];
    pl.drawEllipse(pl.point(p), 2.0, 2.0);
  }

  pen.setColor(Qt::red);
  pen.setStyle(Qt::DashLine);
  pl.setPen(pen);
  n = m_vertexIdx.size();
  for (int i = 0; i < n; i++) {
    vec& p = m_point[m_vertexIdx[i]];
    vec& q = m_point[m_vertexIdx[(i+1) % n]];
    pl.drawLine(pl.point(p), pl.point(q));
    pl.drawEllipse(pl.point(p), 3.0, 3.0);
  }
}
